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Abstract 

We present preliminary results from multidimensional numerical studies of pair instability supernova (PSN), study- 
ing the fluid instabilities that occur in multiple spatial dimensions. We use the new radiation-hydrodynamics code, 
CASTRO, and introduce a new mapping procedure that defines the initial conditions for the multidimensional runs in 
such a way that conservation of physical quantities is guaranteed at any level of resolution. 
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1. Introduction 

The first stars that formed after the big bang may 
have a characteristic mass scale of around hundred so- 
lar masses (M Q ) HI 0, . The stars with initial mass 
between 140 M and 260 M end their lives in a very 
powerful explosion, called a pair-instability supernova 
(PSN) 0, B 0]. During the evolution of those mas- 
sive stars, after central carbon burning, the core of the 
star reaches a sufficiently high temperature that electron 
and positron pairs can be produced. At this time, ra- 
diation energy turns into rest mass for these pairs. It 
softens the adiabatic index y of the gas below the crit- 
ical value of 4/3 and triggers a rapid contraction that 
leads to explosive burning of oxygen and silicon. The 
energy released raises the pressure enough to turn the 
contraction around into a energetic thermonuclear ex- 
plosion (~ 3 - lOOxlO 51 erg). These supernovae may 
play an important role in the synthesis of heavy ele- 
ments \M and their energetic feedback to their sur- 
roundings can affect the structure and evolution of the 
early universe BEH. Although several PSN candidates 
have been observed ifTol 11], there are still many dis- 
crepancies between models and observations [ 12]. The 
current theoretical models for PSN are all based on one- 
dimensional calculations |@, Q] ; until now, no multidi- 
mensional simulations have scarce. Here we study how 
multi-dimensional fluid instabilities affect the mixing of 
elements and possibly even the overall nucleosynthesis 
and energetics of PSN. 
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In this paper, we first introduce our numerical ap- 
proach in §|2] Then, in §|3] we introduce a new mapping 
procedure that defines the initial conditions for the mul- 
tidimensional simulations in such a way that conserva- 
tion of physical quantities is guaranteed at any resolu- 
tion. We discuss the results of our simulation in §|4]and 
present our conclusions in § 



2. Numerical Approach 

We start our simulations using one-dimensional mod- 
els obtained from the KEPLER code fTHl . spherically 
symmetric Lagrangian code that followed the evolu- 
tion of a 150 M star up to ten seconds before max- 
imum compression. Then we map the resulting one- 
dimensional profiles into 2D and 3D to serve as the ini- 
tial conditions for CASTRO Q. We then evolve the sim- 
ulation for about one hundred seconds. This is period 
during which thermonuclear burning releases almost all 
of the energy of the explosion. 

CASTRO [14] is a new, massively parallel, multi- 
dimensional Eulerian AMR radiation-hydrodynamics 
code for astrophysical applications. Time integration 
of the hydrodynamics equations is based on a higher- 
order, unsplit Godunov scheme. Block-structured adap- 
tive mesh refinement (AMR) and sub-cycling in time 
enable the use of high spatial resolution where it is most 
needed. We refer the readers to 111411511 for details. 

CASTRO allows the user to specify their equation of 
state, reaction network, and method for calculating self- 
gravity. For the simulations presented here, we use the 
Helmholtz equation of state lldl . Our reaction network 
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contains 19 isotopes and 86 reaction rates [13]. This 
network is sufficient for us to model the burning pro- 
cesses and their energy release to sufficient accuracy. 
We use the monopole approximation for self-gravity. 

3. Mapping of Initial Model in Multi-D 

Since our initial model will be in a state close to hy- 
drostatic equilibrium, we need to be careful how we 
map the ID spherically symmetric data given on a non- 
uniform Lagrangian grid to a multidimensional Eulerian 
grid. Here we present a method that numerically con- 
serves quantities such as mass and energy that are ana- 
lytically conserved in the evolution equations. Whereas 
this does not guarantee that the initial data will be in per- 
fect numerical hydrostatic balance, it is at least a physi- 
cally motivated constraint and is sufficient for our simu- 
lations. The algorithm as described below is somewhat 
specific to our data set, but can be easily generalized to 
other cases of mapping ID data to higher dimensions. 

Our original ID data is given as cells of known size 
(radius and mass) with known average values of inten- 
sive quantities (density, internal energy). Velocity is 
only given for the zone boundaries and hence total mo- 
mentum conservation is somewhat arbitrarjfl as is ki- 
netic energy conservation. 

The first step is to construct a continuous function 
(C°) that conserves the physically conserved quantities. 
We found an ideal choice is to use "volume coordinate", 
V, the volume enclosed by a given radius from the cen- 
ter of the star. We then use densities (mass density, en- 
ergy density) such that the integral under the curve 
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is the total amount of the quantity X with space density 
px between the volume coordinates V\ and V%. 

Here we use a piecewise linear function that preserves 
the original monotonicity of the data (does not create ar- 
tificial extrema) and is bounded by the original extrema 
of the data (Fig. [TJ. The scheme causes some "smear- 
ing" (smoothing) of the data which, however, is limited 
to less than one zone width of the original data. We use 
a two-step process: First we construct a linear interpo- 
lation across the interface between two zones which ex- 
tends to the half-width of the smaller of the two zones. 
What is cut off from one zone by the interpolation func- 
tion (a or b) is added to the interpolation in the neighbor- 
ing zone (a' - a and V = b). This would usually result 




1 since the initial model is spherically symmetric its total physical 
momentum is always zero for every radius bin 



Figure 1: Schematic for constructing a conservative density profile. 



in two "kinks" (change in slope) inside a zone (middle 
zone in Fig. [TJ with a flat piece that usually is a poor 
approximation to the average gradient. We now correct 
the interpolation within each zone such that there is only 
one "kink" by finding a point in the flat piece such that 
we have equal area enclosed by the triangles on either 
side connecting from the new point to the value of the 
previous interpolation function at the boundaries of the 
zone (c = c')- 

In principle one could now integrate this function 
over the volume of the cell of the target grid to obtain the 
total amount of X within this cell. For the sake of gen- 
erality of the interpolation function, we use an adaptive 
iterative sub-sampling method to obtain an acceptably 
converged integral: We evaluate the interpolation func- 
tion at one point and multiply by the zone volume, then 
subdivide the zone and sample the function for each 
of these sub-volumes, multiply by the sub-volumes and 
sum up the result; this is repeated recursively for each 
sub-zone until the results from the last two levels agree 
to within the desired relative accuracy. Table 1 shows 
the results of mass deviation from original data after 
initial mapping by using linear interpolation and our al- 
gorithm; the number inside the () indicates the resolu- 
tion of the mapped data. By comparing both results, our 
method has much higher accuracy and less depends on 
resolution and spacial dimensions. 

This method is used to map internal energy density 
and the partial densities of the chemical species. The 
total density (or mass within the zone) is then obtained 
from the sum of partial densities (masses); pressure and 
temperature are obtained from the equation of state. 
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Table 1 



Dimension Linear Interpolation Our Algorithm 

ID "2.91%(1024) ~a000254%(1024) 

2D " 47.45%(1024 z ) ~~0l)05%(1024 2 ) 
3D 127.05%(128 3 ) 1.85%(2 3 ) 



4. Results and Discussion 

We performed our simulations on Calhoun at the 
Minnesota Supercomputing Institute (MSI). The 2D 
run used about 4, OOOCPUh; the 3D run used over 
50, 000 CPUh. The results presented here are from 2D 
runs of 150 M Q PSN in cylindrical symmetry where we 
simulated only one hemisphere. 

Figure [2^ shows the specific nuclear energy gener- 
ation rate in the inner 2x10 10 cm domain at about 60 
sees after maximum compression. The color coding 
shows the specific nuclear energy generation rate and 
in units of ergs/s/g on a logarithmic scale. We find 
Rayleigh-Taylor (RT) instabilities develop at the edge of 
the oxygen-burning shell. Later these instabilities will 
grow further an affect such properties as the observable 
light curve of the supernova. Figure [2J) is a close-up of 
the RT instability. 

5. Conclusions 

We have presented preliminary results from our first 
multidimensional numerical study of the evolution of 
pair instability supernova using the new Eulerian AMR 
radiation-hydrodynamics code CASTRO. We simulated 
the formation of Rayleigh-Taylor instabilities in the ex- 
plosion of a 150 M pair-instability supernova. We have 
introduced a new mapping method that can be used to 
define the initial conditions for multidimensional simu- 
lations from one-dimensional initial data in such a way 
that conservation of physical quantities, monotonicity, 
and continuity are guaranteed at any resolution. 
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Figure 2: Energy generation rate 
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